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ABSTRACT 

The generation and propagation of small amplitude 
acoustic waves in a homogeneous, loss-free, compressible 
fluid is studied by the finite element method. A diaphragm 
mounted in an infinite rigid baffle generates acoustic 
waves in a semi-infinite fluid region. Steady-state pres- 
sure distribution is found for a hemispherical region with 
boundary reflection suppressed through use of a radiation 
condition. The computer program developed for the purpose 
utilizes iso-parametric finite elements with curvilinear 
boundaries. Incorporated in the program is a versatile 
mesh generator which minimizes the quantity of input data. 
Acceptable agreement with analytic results is obtained when 
there are at least four elements per wave-length. 
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I. INTRODUCTION 



The field of acoustics plays an important role in the 
development of the complex sonar systems used in the United 
States Navy today. It is the intent of this paper to open 
a new avenue to the study of the active sonar system. 
Attention is confined to the transducer and the periodic 
pressure field induced in a fluid medium by harmonic 
oscillations of the transducer piston or diaphragm. 

The study is restricted to the steady-state performance 
of a transducer mounted in an infinite rigid baffle and 
radiating into a half-space filled with a homogeneous fluid. 
Both the pressure field and the acoustic impedance of the 
transducer are determined. Figure 1 illustrates the 
geometry of the region. 




Figure 1. Geometry of Region. 
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The method used in this study is the two-dimensional 
finite element technique. Since it is manifestly impossible 
to represent an infinite region with a finite number of 
finite elements, the stratagem of the anechoic chamber is 
simulated mathematically. This is accomplished by using 
radiation (non-reflection) boundary conditions. 

The diaphragm, being circular, forms a surface of 
revolution when deflected. Thus, the instantaneous pres- 
sure distribution has rotational symmetry about an axis 
through the diaphragm center and normal to the rigid baffle 
(z axis of Figure 1). A finite element assembly representing 
the region of interest may be constructed by dividing the 
corresponding portion of the r~z plane of Figure 1 into sub- 
areas as shown in Figure 2. 




Figure 2. Finite Element Grid 
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In Figure 2, an individual rectangular element represents 
the cross-section of a toroidal element generated by revolu- 
tion about the ?■ axis. The rotational symmetry permits 
calculation of element properties using two-dimensional 
technique by simply taking the thickness (normal to the 
plane of the figure) proportional to r. 
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II . F INITE ELEMENT FORMULATION OF 
THE GOVERNING EQUATION 



A. EQUATION OF FLUID MOTION 

The governing partial differential equation for the pres- 
sure p in a homogeneous, inviscid, compressible fluid as 
used in acoustics (e.g., see Ref. [1]') is 

v 2 P = p (1) 

c 

2 

where V is the Laplacian operator, c is the acoustic 
velocity, and superior dots represent differentiation with 
respect to time. For the purpose of the present study, 
boundary conditions take the form 



9jp 

8n 




( 2 ) 



where the coordinate n is measured along the normal to the 
boundary surface, p is the fluid density, and v is the 
component of fluid particle velocity in the n direction. 

B. SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS 

The finite element method replaces Eqs. (1) and (2) by 
a set of simultaneous ordinary differential equations with 
independent variable time and dependent variables consisting 
of pressures at a large number of discrete points (nodes) 
of the region. The region is divided into a number of 



^Numbers in brackets designate references listed on page 68. 
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sub-regions called "elements." Within each element it is 
postulated that the pressure p may be expressed in terms of 
the pressure p^ at the n nodes of the element. This rela- 
tion is expressed as 

n 

P = 2 N.p. (3) 

i=l 

where the "shape" (or interpolation) functions N. are 
functions of the spatial coordinates. If the nodes are 
located at suitably chosen points on the element surface, 
unique sets of shape functions will guarantee the necessary 
continuity of pressure p along interelement boundaries as 
pointed out by Zienklewicz [2]. The elements used in this 
study are two-dimensional eight-noded iso-parametric 
elements. Further details are given in Appendix C. 

It has been shown by Zienkiewicz and Newton [3] that 
the finite element equation for k nodes, m of which are 
located on the fluid-diaphragm interface, may be written as 

t Q 3 { p } + [D] {p } + [ H] (p } = - p[L] T {6) (H) 

where {p} is a k xl vector of nodal pressures, [Q] , [D] , and 

[H] are kxk symmetric matrices, {6} is an m x 1 vector of 

T 

diaphragm displacements, and [L] is a k xm matrix. Elements 

T 

of [Q] , [D] , [H] , and [L] are all real constants. Formulas 

for calculating these matrices are given in Appendix E. The 

. T •• 

terms [D] (p) and [L] {6} represent boundary conditions and 

will be discussed in following sections. 
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C. DIAPHRAGM BOUNDARY CONDITION 



The diaphragm displacement 6 is a prescribed function of 

time. It can be represented by choosing a number of points 

(nodes) on the diaphragm and using a set of shape functions 

n! to interpolate the local displacement as 
<3 

m 

6=1 N.6 . (5) 

J-l J J 

where the 6 . are the displacements at the m diaphragm nodes. 

J 

The boundary condition of Eq . (2) is formulated in terms 

of the normal component of fluid particle velocity. At 
the diaphragm surface this must, in the absence of cavita- 
tion, be equal to 6. Using this equality, the fluid- 
diaphragm interface condition gives the term 

- P [ L ] T { 6 > 

in the finite element equation (Eq. ( ; l)). The resulting 
vector has nonzero elements only for those pressure nodes 
located on the diaphragm surface. It is noted that the 
boundary condition on the rigid baffle is a special case, 
in that the particle velocity = 0 , which implies v fi = 0 . 

D. RADIATION (NON-REFLECTION) BOUNDARY CONDITION 

A normally incident wave meeting the interface bet ween 
two media will not be reflected if the specific acoustic 
impedances of the two media are equal. The specific acoustic 
impedance z is defined to be the ratio of pressure to 
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particle velocity. For plane harmonic waves this is a real 
quantity and the pressure can be given as 



p = pcv 

i k n 



( 6 ) 



for a wave progressing in the positive n direction. Differ- 
entiating F.q . (6) with respect to time and using this result 

to eliminate v^ from the boundary condition (Eq. (2)) gives 



2R = 

9n 




(7) 



Equation (7) then becomes the boundary condition for non- 
reflection of normally incident plane waves. The term 

[D] {p} 

in Eq . ( *1) results from this criterion. This vector has 

nonzero elements only for those nodes which lie on the non- 
reflecting portion of the boundary. 
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III. GEOMETRY OP THE REGION 



The initial solution of the problem using the geometric 
configuration of Figure 2 resulted in significant error in 
the pressure distribution when the solution was checked 
against known results. The principal cause of this error 
was believed to be the shape of the non-reflection boundary. 
By the nature of the geometry, the harmonic oscillating 
diaphragm will create progressive pressure waves which 
propagate outward and strike the boundary of Figure 2 at 
various angles of incidence. The particle velocity, v , 
will not always be normal to the non-reflection boundary 
and the boundary condition required by Eq. (2) will not be 
satisfied. Thus, a modification to the representation of 
the boundary to provide substantially normal incidence led 
to the geometry illustrated in Figure 3. For this modified 
radiation boundary , in three dimensions the surface is 
hemispherical . 

In order to represent adequately the circular arc of 
the non-reflection boundary, two-dimensional iso-parametric 
finite elements are used. The element thickness normal to 
the plane of the figure is directly proportional to the 
radius r. Iso-parametric elements are further discussed 
in Appendix C, 

Note the distortion of the grid in Figure 3- The reason 
for this distortion and the method of generating the grid 
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r 




Figure 3* Modified Finite Element Region 



pattern is discussed in Appendix D. The distorted 
"rectangles" still represent the cross sections of toroidal 
elements . 
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IV. WAVE THEORY AT RADIATION BOUNDARY 



A one-dimensional study of a pulsating sphere acting on 
a fluid medium resulted in a more accurate solution when 
spherical wave theory at the non-reflection boundary was 
introduced. Based upon the success of these results, a 
further modification was made to the computer program 
(Appendix B) to incorporate spherical wave theory in the 
two-dimensional study. This was accomplished by developing 
the non-reflection boundary condition for spherically diverg- 
ing waves as shown below . 

The specific acoustic impedance for spherically diverg- 
ing waves is given in Ref. [1] as 



Z 



p 

v 

n 



= pc 



to 

— r 
c 

(0 

c r '-- J 



( 8 ) 



where to is the circular frequency and r is the radial 
distance from the "center." The particle acceleration v n 
can be written as 



v 

n 



J MV 



n 



( 9 ) 



which allows the boundary condition to be represented by 



8r 



Pi MV 

Kd n 



( 10 ) 



From Eqs . (8) and (10), the boundary condition is 



( 11 ) 



8n 

3r 



,i mp 
c 





or 




( 12 ) 



The imaginary term on the right hand side of Eq . (12) is 

represented by the term [D]{p} = jw[D]{p) in Eq . ( 4| ) . Thus, 

the contribution of the real term can be written as 



and can be added to the system of equations given in Eq . ( J 4). 

The spherical contribution is only added to the nodes 
which lie on the non-reflection boundary, thus providing a 
spherical wave approximation to the plane wave theory at 
the non-reflection boundary. 



§[DHp} 



(13) 
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V . COMPUTER SOLUTION TECHNIQUE 



Assume the diaphragm is oscillating with a harmonic dis- 
placement represented in the complex form by 

{6} = {S o }e jaJt (14) 

where { 6 q } is a vector of prescribed displacement amplitudes 

at the m nodes of the diaphragm-fluid interface. Elements . 

of (6 } are real. The pressure in the fluid medium is then 
o 

represented by 

(P) = (P o )e Jwt (15) 

where {P q } is a kx 1 vector of pressure amplitudes (complex) 
at the k nodes of the region. 

Inserting Eqs . (14) and (1| ) into Eq . (4) gives 

([H] - w 2 [Q] + j w [ D ] ) { p } = w 2 {R) (16) 

where (B) - [L] T {6 q }. 

In the computer program, the system matrices [H] , [OJ , 

[D] , and (B) are assembled from the contributions of each 
individual finite element in various subroutines. Since 
only a small number of nodes are located on the non- 
reflection boundary, the matrix [D] is sparse. In order to 
reduce computer storage requirements, [D] , which is symmetric, 
is stored as three vectors: the principal diagonal and the 

two diagonals immediately above. At the same time, the 
approximated spherical contribution j-[D] at the 
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non-reflection boundary is computed. These quantities are 
real and are added directly to the matrix [II] in the appro- 
priate locations. 

The final subroutine introduces the value of omega and 
multiplies the matrices [Q] , [D] , and {B} by the appropriate 
power of omega before the left-hand side is summed to form 
one complex matrix. Thus the system of equations is 
reduced to 

[S]{p] = {F} (17) 

where [S] is a k xk complex matrix and {F} is a column 
vector of order k x 1 in which only m elements are nonzero. 

The complex pressure amplitudes in Eq . (17) are found 
by using the Gaussian elimination and back substitution 
method . 

The absolute value of the pressure and the phase angle 
are calculated in the program and printed out in array 
form. Additionally , a system energy balance and the 
diaphragm acoustic impedance are determined. 

The program provides for multiple-problem solution by 
using incremented values of omega so that a range of 
frequencies can be studied. 
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RESULTS 



VI . 

Problems studied were chosen to verify the accuracy of 
the method through comparison with available, exact 
(analytic) solutions. Results for each problem are 
discussed below. 

Problem 1, Piston in a Rigid Tube 

To verify the correctness of the computer program, 
finite element solutions were obtained for the problem of 
a harmonically oscillating circular piston at one end of 
a fluid-filled, rigid pipe of equal diameter. An infinite 
length was simulated by terminating the finite element 
region with a plane non-reflection boundary at a distance 
of one piston diameter. This is, of course, a one dimen- 
sional problem with uniform pressure amplitude throughout 
and phase angle is a linear function of distance along 
the pipe . 

Detailed results are not presented for this essentially 
trivial test problem. It was found, however, that there 
was excellent agreement with theory for long wave-lengths . 
For wave-lengths less than approximately four times the 
element width (measured parallel to the direction of wave 
propagation), the accuracy deteriorated. This is under- 
standable, since the shape functions employed represent 
the instantaneous pressure profile across an element by a 
parabolic arc. The lower limit on wave-length corresponds, 
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for a given element grid, to an upper limit on frequency. 
Problem 2 . _ Piston in an Infinite Baffle 

The problem explored most extensively is that of a 
harmonically oscillating circular piston surrounded by an 
infinite rigid baffle and radiating into a semi-infinite 
fluid region. For this problem, the ra.diation boundary is 
a hemispherical surface with the center located at the 
piston center and the radius is four times that of the 
piston. The element grid. (see Figure Appendix A) contains 
16 elements. It is derived by mapping a x square array 
from the plane into a quadrant of a circle. 

Analytic results for this problem are readily avail- 
able (e.g., see [1], [^)). In Table I, the pressure ampli- 
tudes and phase angles along the z axis, obtained by the 
finite element method (FEM), are compared with exact 
results. Tw o sets of finite element results are given. 

For those designated FEM1 , the spherical wave correction at 
the radiation boundary is included. For FEM2, the correc- 
tion is omitted. Details are given in Table I for two 
circular frequencies, w = .1 and .6 rad/sec. The corres- 
ponding wave-lengths are approximately 63 ft and 10 feet. 

Examination of Table I shows that FEM1 results give 
generally good agreement with exact results. Those for 
FEM2 are considerably poorer, although the discrepancies 
are smaller for w = .6 rad/sec. This is understandable since 
the 1/r contribution to specific acoustic impedance is 
relatively less important as w increases (see Eq . (12)). 



19 



It is believed that a more refined finite element mesh 
(more elements) should afford significant improvement. 

Since the computer program for this problem required 
320K bytes of core storage, it was not feasible to 
investigate this without major program revisions. 

In predicting transducer performance the resultant fluid 
reaction on the piston is needed. This is normally reported 
in the form of piston radiation impedance, which is complex 
and is found by dividing the piston force by the piston 
velocity. The real part of the radiation impedance is 
called the radiation resistance and the imaginary part is 
called the radiation reactance. 

Table II provides a comparison with exact values for 
the radiation resistance and reactance. Results are given 
for FEM1 and FEF2 over a range of frequencies. Again it is 
observed that the spherical correction significantly improves 
the results. Note that the radiation resistance obtained 
by FEM1 continues to give good results through to = .9 rad/ 
sec , 

Computer solutions of this problem were obtained for a 
range of frequency from to = .1 rad/sec to to = 1.2 rad/sec 
by .1 rad/sec increments. Comparisons at the upper end of 
this range are omitted from Tables I and II because the 
departure from exact results does not permit useful employ- 
ment of the finite element method without grid refinement. 
Complete computer results for pressure amplitude, phase 
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TABLE I 

PRESSURE AMPLITUDE AND PHASE ANGLE ON AXIS 







PRESSURE AMPLITUDE 


PHASE ANGLE 


(JO 


z 


EXACT 


FEM 1 


FEM2 


EXACT 


FEM1 


FEM2 


rad/sec 


ft 


lb/ft 2 


lb/ft 2 


lb/ft 2 


deg 


deg 


deg 


. 1 


0 


.01997 


. 02061 


. 01794 


-5.7 


-5.7 


- 11.6 


. 1 


2 


. 00828 


, 00889 


.00676 


- 13.8 


-13.1 


-31.9 


.1 


4 


. 00472 


.00469 


.00376 


-24 . 3 


-24.9 


-68.3 


.1 


6 


. 00325 


.00329 


.00337 


-35.3 


-35.8 


- 92.1 


. 1 


8 


. 00246 


. 00247 


. 0033^ 


-46 . 5 


-47.4 


-107.3 


.6 


0 


. 678 


.699 


. 625 


-34 . 4 


-33.2 


-34.7 


.6 


2 


. 295 


. 300 


.273 


- 83 .O 


-78.7 


-8o . 2 


. 6 


4 


.169 


. 172 


.186 


-145.6 


-145.4 


- 156.5 


.6 


6 


• 117 


. 116 


. 106 


-211.8 


-211.0 


- 210.6 


. 6 


8 


.089 


.086 


.084 


-279.3 


-279 • 2 


- 289.3 



TABLE II 

PISTON RADIATION RESISTANCE AND REACTANCE 



w 

rad/sec 


RESISTANCE 


REACTANCE 


EXACT 
lb . sec 
ft 


FEM 1 
lb . sec 
ft 


FEM 2 
lb . sec 
ft 


EXACT 
lb . sec 
ft 


FEM1 
lb . sec 
ft 


FEM 2 
lb . sec 
ft 


. 1 


.249 


. 254 


.451 


2 . 11 


2.06 


1.68 


.2 


.979 


• 978 


. 687 


4 . 09 


3.97 


3.56 


.3 


2.13 


2.12 


1.43 


5.81 


5.65 


5.83 


.4 


3.61 


3-62 


3-75 


7.18 


6.99 


8.11 


.5 


5.32 


5.34 


6.75 


8.13 


7.87 


7.43 


.6 


7 . 12 


7.12 


6.6 1 


8.62 


8.31 


7.00 


.7 


8.89 


8.87 


7.50 


8.67 


8.35 


8.72 


.8 


10. 51 


10.53 


10.89 


8.32 


7.99 


9.59 


• 9 


n. 90 


11.97 


13.60 


7.64 


7.24 


6 . 69 



FEM1 - Finite element results with spherical correction. 
FEM2 - Finite element results without spherical correction. 

Phase angle given is for pressure relative to piston 
displacement . 

Fluid density = 1 lb.sec/ft^ Piston Radius = 2 ft . 

Acoustic Velocity = 1 ft/sec. Piston Amplitude = 1 f t . 
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angle, and impedance components at several frequencies are 
given in Appendix A. 

Problem 3- Parabolic Diaphragm in an Infinite Baffl e 

This problem, a modification of the previous problem, 
replaces the piston by a circular diaphragm having a fixed 
edge and having a deflection surface in the form of a 
paraboloid of revolution. Study of this problem was under- 
taken in an effort to improve agreement in the energy 
balance check (discussed below) incorporated into the 
program. It was believed that elimination of fluid 
velocity discontinuity at the edge of the piston in problem 
1 might result in a significant improvement. 

Analytic results for this problem are not available in 
standard acoustic texts, but conventional techniques can 
be utilized to obtain the pressure p along the z axis in 
complex form as 



P 



_ 2j pc^o 



ua 




2 2 \ 

, , u a . wz \ 

i + — 2 r 

2 c C / 



. wz 



- 1 - 



J tot 



( 18 ) 



where 



S 0 = displacement amplitude at diaphragm center 
a = diaphragm diameter 
s = "slant" distance = ( z^ + a^)^ 



For this problem only a single finite element solution 
(FEM3) was obtained. The spherical correction is included 
Table III gives a comparison with exact results (Eq. (18)) 
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TABLE III 



PRESSURE AMPLITUDE AND PHASE ANGLE ON AXIS 



w 

rad/sec 


3 N 

i 


PRESSURE AMPLITUDE 


PHASE 


ANGLE 


EXACT 

lb/ft 2 


FEM3 

lb/ft 2 


EXACT 

deg 


FEM3 

deg 


.1 


0 


.01322 


. 01*127 


-*l • 3 


-*4 . 1 


. 1 


2 


.00*438 


. 00*l6l 


-13.1 


-12.6 


.1 


*1 


. 002*1 0 


.002*12 


-23-7 


-2*1 .2 


.1 


6 


.0016*] 


.00166 


-3*1 . 8 


-35.5 


.1 


8 


.0012*] 


.00125 


— *46.2 


-*17.1 


. 6 


0 


. A ] 59 8 


. 493 1 * 


-25.7 


OO 

on 

C\j 

1 


.6 


2 


.1565 


.1590 


-78.5 


- 75.5 


.6 


*1 


.086*] 


. 0867 


-1*12.97 


-1*12.3 


.6 


6 


.0589 


.058*1 


-209 • *4 


-209 . *1 


.6 


8 


.0*1*] 5 


.0*136 


-277.8 


-278.1 



FEM3 - Finite element results with spherical correction. 

Phase angle is for pressure relative to diaphragm displacement. 

2i 

Fluid density = 1 lb.sec/ft . 

Acoustic velocity = 1 ft/sec. 

Diaphragm radius = 2 f t . 

Diaphragm amplitude (center) = 1 ft. 



for pressure amplitudes and phase angles for values of 
w = .1 and .6 rad/sec. 

Examination of Table III shows that agreement near the 
diaphragm (small values of z) is not as good as FEM1 in 
Table I. At the radiation boundary however, the FEM3 errors 
are comparable to those for FEM1. 
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Energy Balance Check 



Prior experience, by others, suggested that an energy 
balance check would be a useful feature of a finite element 
program. Specifically, since solutions represent steady- 
state behavior, the input power, averaged over a cycle, 
should equal the corresponding average value of the output 
power. Input occurs at the piston, or diaphragm, and output 
at the radiation boundary. Expressions for these quantities 
are derived and incorporated in the computer program. 

A number of program errors in the energy balance calcu- 
lations (Appendix B) prevented valid results until the 
final stage of the investigation. Thus, the values of 
average power in and average power out are not included in 
the computer output (Appendix A). Recent computations, 
with program errors corrected, have shown excellent agree- 
ment of these quantities to at least six significant digits 
in the lower range of the circular frequency. Even at 
higher frequencies , where solution accuracy (pressure 
amplitudes and phase angles) deteriorates markedly, five 
significant digit agreement is maintained. 



24 



VII . CONCLUSIONS 



It has been shown that the finite element method can 
provide acceptable results concerning generation and prop- 
agation of acoustic waves. The parabolic iso-parametric 
element provides an economical, accurate means of repre- 
senting regions having curved boundaries. 

The principal limitation of the method is evident from 
the empirically deduced results that, for acceptable 
accuracy, a single element must not span more than one- 
quarter wave-length . For wave-lengths of the order of 
magnitude of the piston radius, the number of nodes 
required would be around 900. For the present program, 
which already utilizes 320K bytes of core storage for 65 
nodes, this is far beyond capability. 
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APPENDIX A 



COMPUTER OUTPUT 

The rows of pressures and phase angles shown in the 
computer output correspond to the numbering system shown on 
the finite element region. Figure 

The zeros printed in the array of pressures and phase 
angles do not indicate computed values, but facilitate 
printing the output in array form. 



Row 

Row 

Row 

Row 

Row 




Figure Format for Pressure and Phase Angle Array 



26 



The output data are arranged according to the following 
sequence . 

FEM1. Piston in infinite baffle with spherical correction 
added at the non-reflection boundary. 

FEM2 . Piston in infinite baffle without spherical correc- 
tion added at the non-reflection boundary. 

FEM3- Diaphragm in infinite baffle with spherical correc- 
tion added at the non-reflection boundary. 
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COMPUTER OUTPUT-GENERAL 
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IF(I.EQ.l) K=1 

CALL DRAW(NNX,XDRAW, YDRAW,K»0» LABE L » I T I TLE t 0 , 0, 0 , 0 » 0 * 0 , 5 , 7 , 0 
I=I+(4*NXDIV+2 > 
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IF(C0ABS( A( I ♦ J ) ) -EQ. 0.D0) GO TO 404 
A( I,J)=A( I 1 J ) / A ( I, I> 

DO 406 K= J j N 

A< K » J ) =A( K, J I -A( K» I ) *A ( I , J > 
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APPENDIX C 



ISO-PARAMETRIC ELEMENTS 

The family of iso-parametric elements, recently developed 
by Zienkiewicz , Irons, and others, has brought new versa- 
tility and power to the finite element method. The special 
usefulness of these elements is that they may be shaped to 
conform to curvilinear boundaries. Attention is confined 
here to txvo-dimensional parabolic elements. A full descrip- 
tion of the entire family is given in Ref. [2] . 

Consider the rectangular element of Figure 5. The 
element has four corner nodes and four midside nodes, 
numbered as shown. Points in the interior, or on the 
boundary, may be located in the local (£,n) coordinate 
system. The dependent variable, pressure p in the present 
application, is defined at any point by Eq . (3), repeated 

here 

8 

P = E N p (3) 

i = l 1 1 

where 

N ± = }(i+ nn i )(i + CC i )(nn i + -D 

for corner node i, and 

N ± = |(l- n ? K ± 2 - C 2 n i 2 )(l + nn. 
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4 <K=i 



n=~l 



Figure 5. Element Configuration in Local Coordinates 



for mid-side node i, 

( ^i j ) are coordinates of node i, 

and p^ is the pressure at node i. 

Note that each shape function has the value unity at node 
i , and is zero at all other nodes . 

A curvilinear element, such as that of Figure 6, may 
be produced by mapping the rectangle of Figure 5 into the 
x,y plane by the relations 

8 8 

x = E N. x. , y = E N.y. (19) 

i=l i=l 

where (x^,y^) are the cartesian coordinates of node i. 

The shape of the element in the x,y plane is uniquely 
determined by specifying the coordinates (x^,y^) for the 
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y 




Figure 6. Curvilinear Elements Mapped 
in Cartesian Coordinates. 

eight nodes. The edges of this curvilinear element are 
parabolic . 

Adjacent elements ’’share" the same three nodes at their 
interface. When mapped to cartesian coordinates, these two 
adjacent elements then have the same curvilinear shape and 
the same x and y coordinates at nodes. This provides con- 
tinuity in the pressure distribution across element 
boundaries . 

The region developed for the present study necessitates 
the use of iso-parametric elements to provide the circular 
non-reflection boundary. This boundary is constructed 
using parabolic arc segments which are mapped by use of the 
shape functions to "fit" a quarter circle boundary m car- 
tesian coordinates. 
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APPENDIX D 



FINITE ELEMENT MESH GENERATOR 



Two subroutines, included in the computer program of 
Appendix B, should prove useful in a variety of finite 
element problems. Together, these subroutines subdivide 
the region, systematically assign element numbers and 
node numbers, calculate nodal coordinates, and produce 
a computer plot of the region to be studied. As an example 
consider the region of Figure 7. 

The region is taken to be one large iso-parametric 
element (a mapping of the rectangle of Figure 5)« Its shap 
is prescribed by specifying the (x,y) coordinates of the 




Figure 7. Basic Geometry of Finite Element Region. 



eight nodes. Subdivision of the element can then be 
accomplished at predetermined values of £ and n by using 
Equation (19). 

Subroutine NODE (Appendix B) requires two input 
parameters , the number of divisions desired in the x direc- 
tion and the number of divisions desired in the y direction 
These parameters are used to calculate the increment values 
of £ and n for subdividing the element into smaller 
elements. The subroutine starts numbering global nodes at 
the upper left hand corner (£=-1, r, = l) then adds increment 
values of £ and assigns global node numbers from left to 
right on the boundary n=l. Similarly, after incrementing 
r) , the second row of global node numbers is assigned, etc. 
Global node numbers are assigned to all corner nodes, mid- 
side nodes and the center of each subdivided element (9 
nodes per element). The center node provides unique 
identification for the element and simplifies the printing 
of computer output in a readable form. Subsequently all 
nodes are renumbered, omitting the center nodes. 

The cartesian coordinates corresponding to the global 
node numbers are calculated by a similar process in the sub 
routine XYPORM and are used in the program to plot the 
finite element grid. Figure 8 illustrates a computer plot 
of the finite element mesh generated for a quarter circular 
region (Figure 7) with an eight foot radius. The number 
of subdivisions in both the x and y direction is equal to 
four. The plotting routine produces straight segments 



65 



y 




Figure 8. Finite Element Mesh Generated 
by Subroutine XYFORM. 

between nodes. In this instance, the element boundaries 
are all parabolic arcs except for those lying on the x and 
y axes. 
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APPENDIX E 



MATRIX ELEMENT FORMULAS 



Individual elements of the matrices [Q] , [D] , [L] , and 
[H] of the system of ordinary differential equations (Eq. (4)) 
are given by Ref. [3] as 




V 

c Re 



N.N.dR 
i J 



d 



ij 




N.N.dS 
i J 



l 



/ N ! N . dS 
Se 1 J 





3N . 

l 

3x 



3N . 3N. 

__i + _i 
3x 3y 




3N. 

i 

3z 




dR 



where Re and Se are the element region and external boundary 
respectively . 

The shape functions N^ , N^ are specified in terms of 
local coordinates and it is advantageous to perform the 
integrations in these coordinates. Gaussian quadrature is 
employed. Details of the procedure are given in Ref. [2]. 
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